home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / LIBIP / Descript.c.old < prev    next >
Text File  |  1999-09-11  |  9KB  |  364 lines
  1. /* 
  2. ** descript.c
  3. ** 
  4. ** Practical Algorithms for Image Analysis
  5. ** 
  6. ** Copyright (c) 1997 ????
  7. */
  8.  
  9. /*
  10. ** DESCRIPT
  11. **
  12. **  evaluation of Fourier descriptors for a plane curve
  13. **    encoded in a set {delta_phi(k), l(k)};
  14. **
  15. **  references:
  16. **  C. T. Zahn, 
  17. **  Proc. Joint Intern. Conf. on Artif. Intelligence, May 1969, pp. 621 - 628;
  18. **  SLAC Report 70, Stanford Linear Accelerator Center, 1968;
  19. **  C. T. Zahn & R. Z. Roskies, IEEE Trans. Comp., Vol C-21, 269 - 281 (1972);
  20. **
  21. **  -->the truncated Fourier expansion derived by Z & R is computed 
  22. **     in form of inner products using array processor functions
  23. **
  24. **  -->test shapes (see Z & R): 
  25. **        four (xdrawfour.c), inverted L (fdzahnl.c)
  26. **
  27. ** ---------------------------------------------------------------------------
  28. **
  29. */
  30. #include <stdio.h>
  31. #include <stdlib.h>
  32. #include <math.h>
  33. #include "ip.h"
  34.  
  35. #define    EPS        0.00000001
  36. #define    SQ2        1.414213562
  37.  
  38. #define    ON        1
  39. #define    OFF        0
  40. #define    FD_DEBUG    ON
  41.  
  42.  
  43. /*
  44.  * descriptors()
  45.  *   DESCRIPTION:
  46.  *     evaluate ZAHN-ROSKIES Fourier descriptors
  47.  *   ARGUMENTS:
  48.  *     dphik: delta phik for polygon (float *)
  49.  *     dlk: delta lk for polygon (float *)
  50.  *     mcp: length of dphik and dlk (long)
  51.  *     a_n: Fourier coefficients array (float *)
  52.  *     b_n: Fourier coefficients array (float *)
  53.  *     n_order: order of Fourier descriptors to compute (long)
  54.  *   RETURN VALUE:
  55.  *     none
  56.  */
  57.  
  58. void 
  59. descriptors (float *dphik, float *dlk, long mcp, float *a_n, float *b_n, long n_order)
  60. {
  61.     int k, n_ord;
  62.     int fund_flag = OFF;
  63.     
  64.     float real_length, length;
  65.     float dc_term, a, b, part_sum;
  66.     float sne1k, snen, csn1k, csnn;
  67.     float *angle;
  68.     float *sne1, *sne, *csn1, *csn;
  69.     float *mod, *arg;
  70.  
  71.  
  72. /*  
  73. **  compute (in place) running sums of the elements in dlk 
  74. **  and the corresponding phase angles required in the evaluation
  75. **  of the fundamental
  76. */
  77.     angle = (float *)calloc((int)mcp, sizeof(float));
  78.     sne1 = (float *)calloc((int)mcp, sizeof(float));
  79.     csn1 = (float *)calloc((int)mcp, sizeof(float));
  80.     sne = (float *)calloc((int)mcp, sizeof(float));
  81.     csn = (float *)calloc((int)mcp, sizeof(float));
  82.     if((angle == NULL) || (sne == NULL) || (csn1 == NULL) || (sne == NULL)
  83.         || (csn == NULL))
  84.     {
  85.         printf("memory allocation error in function descriptors\n");
  86.         exit(1);
  87.     }
  88.  
  89. /*
  90. ** evaluate contour length, based on curvature points
  91. */
  92.     part_sum = *(dlk+0);
  93.     for(k=1; k<mcp; k++)
  94.     {
  95.         part_sum += *(dlk+k);
  96.         *(dlk+k) = part_sum;        /* array of partial sums */
  97.     }
  98.     length = *(dlk+mcp-1);
  99.         
  100.     for(k=0; k<mcp; k++)
  101.         *(angle+k) = (float)(2*PI*(*(dlk+k))/length);
  102.  
  103.     real_length = (float)(0.5*length);
  104.     printf("\n...length = %f\n", real_length);
  105.  
  106. /*
  107. **  dphik contains elements of delta_phik multiplied by PI/4;
  108. **  to speed up computation, factor out PI 
  109. */
  110.     for(k=0; k<mcp; k++)
  111.         *(dphik+k) *= (float)0.25;
  112.  
  113. #ifdef FD_DEBUG
  114.     printf("...elements of dphik and dlk:\n");
  115.     for(k=0; k<mcp; k++)
  116.         printf("%7.4f  %7.4f\n", *(dphik+k), *(dlk+k) );
  117. #endif
  118.     
  119. /*
  120.  *  compute DC term
  121.  */
  122.     vdot(dphik, dlk, &dc_term, mcp);
  123.     dc_term = (float)(-PI*(1.0 + dc_term/length));
  124.     printf("...and we have a result for the DC term: %f   \n", dc_term);
  125.  
  126. /*
  127. **  execute routine to compute sine and
  128. **  cosine terms for fundamental, then fetch results
  129. */
  130.     for(n_ord=0; n_ord<(int)n_order; n_ord++)
  131.     {
  132.  
  133.     /*  
  134.      *  compute coefficients for fundamental;
  135.      */
  136.         if(n_ord == 0)
  137.         {
  138.             for(k = 0; k < mcp; k++)
  139.                 *(angle + k) = *(angle + k) * (n_ord + 1);
  140.             vsin(angle, sne1, mcp);
  141.             vcos(angle, csn1, mcp);
  142.             for (k=0; k<mcp; k++)
  143.             {
  144.                 csn1k = *(csn1+k);
  145.                 sne1k = *(sne1+k);
  146.                 if( fabs(csn1k) < EPS )  csn1k = (float)0.0;
  147.                 if( fabs(sne1k) < EPS )  sne1k = (float)0.0;
  148.             
  149.                 *(sne+k) = *(sne1+k) = sne1k;
  150.                 *(csn+k) = *(csn1+k) = csn1k;
  151.             }
  152.             fund_flag = ON;
  153.         }
  154.  
  155.  
  156. /*
  157. **  employ double angle formulas to compute higher order coefficients
  158. */
  159.         else
  160.         {
  161.             if( fund_flag != ON )
  162.             {
  163.                 printf("\n");
  164.                 printf("hey, evaluate fundamental first!\n");
  165.                 printf(" something wrong!!\n");
  166.                 exit(1);
  167.             }
  168.             for (k=0; k<mcp; k++)
  169.             {
  170.                 snen = *(sne+k);
  171.                 csnn = *(csn+k);
  172.                 *(sne+k) = snen*(*(csn1+k))+csnn*(*(sne1+k));
  173.                 *(csn+k) = csnn*(*(csn1+k))-snen*(*(sne1+k));
  174.             }
  175. #ifdef FD_DEBUG
  176.             printf("\n...n_ord = %d:\n", n_ord);
  177.             for(k=0; k<10; k++)
  178.                 printf("......sne[%d] = %f   csn[%d] = %f\n",
  179.                         k, *(sne+k), k, *(csn+k) );
  180. #endif
  181.         }
  182.  
  183.  
  184. /*
  185. ** compute coefficients a_n, b_n
  186. */
  187.  
  188.         vdot(dphik, sne, &a, mcp);
  189.         vdot(dphik, csn, &b, mcp);
  190.  
  191.         *(a_n+n_ord) = -a/(float)(n_ord+1); 
  192.         *(b_n+n_ord) = b/(float)(n_ord+1);
  193.     }
  194.  
  195.     printf("...descriptors up to %ld-th order are:\n", n_order);
  196.     for(n_ord=0; n_ord<(int)n_order; n_ord++)
  197.         printf ("%7.3f   %7.3f\n", *(a_n+n_ord), *(b_n+n_ord));
  198.  
  199.  
  200.     mod = (float *)calloc((int)n_order, sizeof(float));
  201.     arg = (float *)calloc((int)n_order, sizeof(float));
  202. /*
  203. **  calculate xy to polar coordinate transformation
  204. */
  205.     vrtop(a_n, b_n, mod, arg, n_order);
  206.     printf("...polar descriptors up to %ld-th order are:\n", n_order);
  207.     for(n_ord=0; n_ord<(int)n_order; n_ord++)
  208.     {
  209.         if( *(arg+n_ord) < 0.0)     *(arg+n_ord) += (float)(2.0*PI);
  210.         printf("%7.3f   %7.3f\n", *(mod+n_ord), *(arg+n_ord));
  211.     }
  212.  
  213.     free(angle);
  214.     free(sne1);
  215.     free(csn1);
  216.     free(sne);
  217.     free(csn);
  218. }
  219.  
  220. /*
  221.  * vdot()
  222.  *   DESCRIPTION:
  223.  *     take the dot product of 2 vectors
  224.  *     and compute the sum
  225.  *   ARGUMENTS:
  226.  *     v1: first vector (float *)
  227.  *     v2: second vector (float *)
  228.  *     vsum: vector sum (float *)
  229.  *     vlen: length of vectors (long)
  230.  *   RETURN VALUE:
  231.  *     none
  232.  */
  233.  
  234. void
  235. vdot(float *v1, float *v2, float *vsum, long vlen)
  236. {
  237.     float rsum = (float)0.0;
  238.     int i;
  239.  
  240.     for(i = 0; i < vlen; i++)
  241.         rsum = rsum + *(v1 + i) * (*(v2 + i));
  242.     *vsum = rsum;
  243. }
  244.  
  245. /*
  246.  * vcos()
  247.  *   DESCRIPTION:
  248.  *     takes the cosine of a vector
  249.  *   ARGUMENTS:
  250.  *     source: input vector (float *) NOTE: radians!!
  251.  *     dest: output vector (float *)
  252.  *     vlen: vector length (long)
  253.  *   RETURN VALUE:
  254.  *     none
  255.  */
  256.  
  257. void
  258. vcos(float *source, float *dest, long vlen)
  259. {
  260.     int i;
  261.  
  262.     for(i = 0; i < vlen; i++)
  263.         *(dest + i) = (float)cos((double)*(source + i));
  264. }
  265.  
  266. /*
  267.  * vsin()
  268.  *   DESCRIPTION:
  269.  *     takes the sine of a vector
  270.  *   ARGUMENTS:
  271.  *     source: input vector (float *) NOTE: radians!!
  272.  *     dest: output vector (float *)
  273.  *     vlen: vector length (long)
  274.  *   RETURN VALUE:
  275.  *     none
  276.  */
  277.  
  278. void
  279. vsin(float *source, float *dest, long vlen)
  280. {
  281.     int i;
  282.  
  283.     for(i = 0; i < vlen; i++)
  284.         *(dest + i) = (float)sin((double)*(source + i));
  285. }
  286.  
  287. /*
  288.  * vrtop()
  289.  *   DESCRIPTION:
  290.  *     converts rectangular cartesian
  291.  *     coordinates to polar
  292.  *   ARGUMENTS:
  293.  *     x: input x vector (float *)
  294.  *     y: input y vector (float *)
  295.  *     mod: modulus output vector (float *)
  296.  *     arg: atan2 output vector (float *)
  297.  *     vlen: vector lengths (long)
  298.  *   RETURN VALUE:
  299.  *     none
  300.  */
  301.  
  302. void
  303. vrtop(float *x, float *y, float *mod, float *arg, long vlen)
  304. {
  305.     int i;
  306.  
  307.     for(i = 0; i < vlen; i++) {
  308.         *(mod + i) = (*(x + i))*(*(x + i)) + (*(y + i))*(*(y + i));
  309.         *(arg + i) = (float)atan2(*(y + i), *(x + i));
  310.     }
  311. }
  312.  
  313. /*
  314.  * msdescriptors()
  315.  *   DESCRIPTION:
  316.  *     wrapper for call to descriptors()
  317.  *   ARGUMENTS:
  318.  *     delta_phik: delta phik for polygon (float *)
  319.  *     delta_lk: delta lk for polygon (float *)
  320.  *     mcp: length of dphik and dlk (long)
  321.  *     a_n: Fourier coefficients array (float *)
  322.  *     b_n: Fourier coefficients array (float *)
  323.  *     n_order: order of Fourier descriptors to compute (long)
  324.  *   RETURN VALUE:
  325.  *     none
  326.  */
  327.  
  328. void
  329. msdescriptors(float *delta_phik, float *delta_lk, long mcp, float *a_n, float *b_n, long n_order)
  330. {
  331.     int k;
  332.  
  333. #ifdef FD_DEBUG
  334.     printf("\nevaluation of Fourier descriptors (Zahn & Roskies)\n");
  335.     printf("--------------------------------------------------\n");
  336.  
  337.     printf("...given: two vectors, delta_phik[] and delta_lk[]\n");
  338.  
  339.     printf(" delta_phik:\n");
  340.     for(k=0; k<(int)mcp; k++)
  341.     {
  342.         printf( "%7.2f  ",*(delta_phik+k) );
  343.         if( (k+1)%8 == 0)
  344.             printf ("\n");
  345.     }   
  346.     printf ("\n   delta_lk:\n");
  347.     for (k=0; k<(int)mcp; k++)
  348.     {
  349.         printf( "%7.2f  ",*(delta_lk+k) );
  350.         if( (k+1)%8 == 0)
  351.             printf ("\n");
  352.     }
  353. #endif
  354.  
  355.     if(n_order > mcp)
  356.     {
  357.         printf("...polygon only has %ld vertices;");
  358.         printf(" cannot yield %ld harmonics!\n", mcp, n_order);
  359.         exit(1);
  360.     }
  361.     descriptors(delta_phik, delta_lk, mcp, a_n, b_n, n_order);
  362.  
  363. }
  364.